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We consider the release from a point source of relatively heavy fluid into a porous sat- 
urated medium above an impermeable slope. We consider the case where the volume of 
the resulting gravity current increases with time like t a and show that for a < 3, at short 
times the current spreads axisymmetrically, with radius r ~ £( a + 1 )/ 4 ) while at long times 
it spreads predominantly downslope. In particular, for long times the downslope position 
of the current scales like t while the current extends a distance W 3 across the slope. For 
a > 3, this situation is reversed with spreading occurring predominantly downslope for 
short times. The governing equations admit similarity solutions whose scaling behaviour 
we determine, with the full similarity form being evaluated by numerical computations 
of the governing partial differential equation. We find that the results of these analyses 
are in good quantitative agreement with a series of laboratory experiments. Finally, we 
discuss the implications of our work for the sequestration of carbon dioxide in aquifers 
with a sloping, impermeable cap. 



1. Introduction 

Horizontal differences in density between two fluids lead to the propagation of so-called 
gravity currents. These currents are of interest in a number of industrial as well as natural 
applications and so obtaining an understanding of the way in which th ey propagate is a 
subject that has motivated a considerable amount of current research ()Huppertl l2006). 

In previous publications, our understandi ng of axisymmetric viscous gravity currents 
on an impermeable b oundary llHuppertlll982|) has been generalised to take account of 
the effect s of a slope jListerHl 992) as well as the p ropagation of a current in a porous 
medium (fHuppert fc Woods! 1 19 9 5: Lvl e et a/.ll2005|) . Here, we consider the propagation 
of a gravity current from a point source in a porous medium at an impermeable sloping 
boundary. Of particular interest is the evolut i on of the current away from the axisym- 
metric similarity solution found by iLvle et al\ ()2005|) . 

We begin by deriving the evolution equations for the shape of a current whose volume 
varies in time like qt a . A scaling analysis of these governing equations reveals the extent 
of the current as a function of time up to a multiplicative constant. The full form of the 
similarity solutions that give rise to these scalings can only be determined by numerical 
means, however, and to do so we modify the numerical code of iListerl ()l992|h For some 
particular values of a, it is possible to make analytical progress; these cases are considered 
separately and provide a useful check of the numerical scheme. We then compare the 
results of the numerical calculations to a series of experiments and find good quantitative 
agreement between the two. Finally, in the last section, we discuss the implications of our 
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(a) 




y = y e (x) 



(b) 




Figure 1. Sketches of a gravity current, of density p + Ap, propagating in a porous medium 
saturated with liquid of density p above an inclined plane, (a) Plan view of the current and (b) 
horizontal section through the current. 

results in geological settings, with particular emphasis on the implications of our work 
for the sequestration of carbon dioxide. 

2. Formulation 



We consider a gravity current consisting of fluid material of density p + Ap in a deep 
porous medium saturated with fluid of density p, which is bounded by an impermeable 
barrier at an angle 9 to the horizontal (see figure for a sketch of the setup). That the 
saturated porous medium is deep in comparison with the vertical extent of the current 
allows us to neglect the motion of the surrounding fluid, simplifying the problem consid- 
erably. We use the natural Cartesian co-ordinate system centred on the mass source and 
aligned with the slope of the impermeable boundary. The depth, h{x, y, t) . of the gra vity 
current is then determined by continuity combined with Darcy's law fsee lBean il988. for 
example) and the assumption that the pressure in the current is hydrostatic, i.e. 



where k is the permeability of the porous medium and p, is the viscosity of the liquid. 
The velocity within the porous medium is therefore given by 



2.1. Governing equations 



P — Pq = ApghcosO — (p + Ap)gz cos9 + pgxsinO (z < h), 
with Pq constant. Here, Darcy's law takes the form 



(2.1) 




(2.2) 




(2.3) 



Using this along with the conservation of mass, we obtain 
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where 4> is the porosity of the porous medium and g' = gAp/p. Equation (|2.4|l is a 
nonlinear advection-diffusion equation for the current thickness, with the two terms on 
the right hand side representing the gravity-driven spreading of the current and its 
advection downslope, respectively. 

It is common to close the system by re quiring that the volume of the current depend 
on tim e like qt a for some constant a > llHuDDertlll98ilListerlll992HHuDDert fc Woods! 
Il995|) . This constraint leads to solutions of self-similar form (as we shall see again in this 
case) but also covers the natural cases of a fixed volume release (a — 0) and a constant 
flux release (a = 1). To impose this volume constraint, (12.411 must be solved along with 

h Ay dx = qt a , (2.5) 

with \y\ = y e (x) giving the edge of the current for x u (t) < x < x n (t). Note that l|2.5|l 
contains an extra multiplicative factor of 4>, which was om itted in the study of an ax- 
isymmetric current in a porous medium bv lLvle et all l)2f)05^ . 

Equations 1|2.4JI and (|2.5(l may be non-dimensionalized by setting T = t/t* , H = h/h* , 
X = x/x* and Y = y/y*, where 

t*=(—p; - — - V " , x*=y* = Vt*, h* =x*iw.0, (2.6) 



4>V 3 tan( 



and 



v ^kpg^n9 
H<f> 

is the natural velocity scale in the problem. In non-dimensional terms, therefore, the 
current satisfies 

along with the volume conservation constraint 

x n f YAX) 

I HdYdX=T a . (2.9) 

X u J-Y e [X) 

2.2. Scalings 

To aid our physical understanding of the spreading of the gravity current, we begin by 
considering the scaling behaviour of the spreading in the limits of short and long times. 
For a < 3, (|2.3|l shows that at short times (T <§; 1) the typical horizontal velocity 
scale is X/T ~ H/X so that H ~ X 2 /T. Further, X ~ Y and volume conservation 
(|2.9|l requir es that HXY ~ T a . From this we therefore find the axisymmetric scalings 
obtained bv lLvle et al\ l)2005|) . namely 

H ~ T 2 ^ 1 , X ~y ^T^r-. (2.10) 

At long times (T> 1), again for a < 3, the typical downslope velocity of the current 
is X/T ~ 1 while in the across-slope direction we have Y/T ~ H/Y. Combined with 
volume conservation XYH ~ T a these scalings lead to 

H~T^, X~T, F~Tf, (2.11) 

so that the current spreads predominantly downslope. It is worth noting here that the 
long time scaling X ~ T is unsurprising because 112.811 may be simpl ified by moving into 
a frame moving at unit speed downslope IjHuppert & Woodslll995IK We also note that 
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Regime Downslope extent Cross-slope extent Thickness 

x y h 



a < 3 t » t* 
a > 3 i«t* 
a > 3 i > t* 



~™ .(^ f - ? ) 1/3 r/3 .(4--) 1/3 ^-3)/3 

, ( Vq \ 1/4 ,(a + l)/4 ( Vq \ 1/4 ,( Q + l)/4 . . ( q tan 9 \ 1/2 , (q- 1) /2 

^tanfly ^0 tan fly ^ y 

Table 1 . Summary of the asymptotic scalings for the dimensions of a gravity current in a porous 
medium at an inclined plane. Here dimensional notation is used for clarity, and t* and V are as 
defined in 12. (il and (12. 71 . respectively. 



the scaling Y ~ 7W 3 j s identical to that found bv lListerl Jl992ft for a viscous current on 
a slope. 

When a > 3, the importance of the two downslope terms (the diffusive and transla- 
tional terms) reverses. In particular, at long times (HHx)x 3> Hx, so that we in fact 
recover the axisymmetric spreading scalings given in 12.100 as being relevant for T>1. 
Conversely, for T< 1 we recover the non-axisymmetric scalings of (|2.11|l . A summary 
of the different scaling regimes expected is given in dimensional terms in tabled 

That we observe axisymmetric spreading if a > 3 and T 3> 1 is surprising, but is 
a consequence of the fact that the downslope flux in a porous medium gravity current 
is only weakly dependent on the local height and so can be swamped by the spreading 
terms in ()2.8(l . In the viscous case, this is not possible because the downslope flux is able 
to remove the incoming flux much more efficiently and penalizes the accumulation of 
material at a particular point more. 



2.3. Numerics 

The axisymmetric sprea ding of a gravity c urrent in a porous medium above an horizontal 
plane was considered bv lLvle et all 1 20051) . In particular, they determined the coefficients 
in the scalings H2.10(l by finding a solution dependent on one similarity variable in this 
case. To determine the prefactors in the non-axisymmetric scaling relations H2.11(l . it is 
necessary to resort to numerical solutions of 12.811 and 12.911 . The numerical code used 
to do this was adapted from that used by iListerl (1992) for a viscous gravity current 
on an inclined plane, with minor alterations to make it applicable to a gravity current 
in a porous medium. This code is an implementation of a finite-difference scheme on 
a rectangular grid with time-stepping performed using an alternating-direction-implicit 
method. Equation l|2.8|l was written in flux-conservativ e form allowing the diffus ive and 
advective terms to be represented by the Il'in s cheme dClauser fc Kiesne"illl987|) . More 
details of the numerical scheme may be found in IListerl l)l992(l . 



3. Special values of a 

In this section, we consider separately particular values of a that are of special interest. 
In some of these cases, it is possible to make progress analytically providing useful checks 
on the numerical scheme discussed in section ROl but they also shed light on situations 
of practical interest. 
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3.1. Constant volume 

As already noted, the differential equation in 12.8fl may be simplified by moving into a 
frame translating at unit speed downslope. However, for general values of a, this cor- 
responds to a point source that is moving uphill in the new frame, complicating the 
analysis. For a current of constant volume, a = 0, there is no distinguished source point 
and we let X' = X — T. The resu lting transformation of (|2.8|) has an axisymmetric 
similarity solution llLvle et alWmi^ . which may be written 

*(W) = ^(^-|^)> (3-D 
where R' = {X' 2 + Y 2 ) 1 ' 2 . 

3.2. Constant flux: A steady state 

For very long times T 2§> 1, we expect that a constant flux current (corresponding to 
a = 1) will approach a steady state, whose shape we now determine. We expect this 
steady shape to be observed far from the nose of the current, since the nose is always 
unsteady, requiring that X <C T. Sufficiently far downstream from the source (X 1), 
the steady shape is given by 

d 2 H 2 _ o 8H 

~dY 2 ~ ~dX 7 {3 - } 

which has a similarity solution of the form H(X, Y) = X^ 1 / 3, f(Y/ X 1 / 3 ) where the 
function / satisfies 

d 2 f 2 If df\ P" 

~'f + V^-)=0, / / dT? = 1, f(± V e) = 0. (3.3) 



drj 2 3 \ drj 
This has solution 



m = l(vt-V% (3.4) 



where r\ e — (9/2) 1 / 3 w 1.651 denotes the position of the current edge in similarity 
variables. 

This shows that far away from the source and nose regions, we should expect the shape 
of unsteady currents to approach Y = (9X/2) 1 / 3 . Superimposing this curve onto the 
numerically calculated current provides a useful check of the numerical scheme described 
in section 1^751 This comparison (see figure [2J shows that, away from both the nose and 
source regions, we do indeed see the steady state shape, though this region is confined 
to T~ x <C X/T <C 1 in the rescaled co-ordinates used in figure [3 

It is interesting to note that the similarity solution (j3.4|l is precisely that given bv lHunpert & Woodsl 
( 1995) for the shape of a two-dimensional current of constant volume spreading in a 
porous medium above an horizontal boundary. This correspondence arises because in the 
steady state case considered here, fluid moves downslope at a constant velocity — inde- 
pendently of its cross-slope position and the current height — so that X is a proxy for 
time. A material slice in the y-z plane thus remains planar as it is advected downslope 
and so spreads laterally in exactly the same way that a fixed volume release does in 
two-dimensions. 



3.3. a = 3 

When a = 3, the non-dimensionalization leading to 1)2. 8JI breaks down because there is no 
longer a characteristic time-scale t* of the motion. Instead, an additional natural velocity 




Figure 2. Numerical evolution of the boundary of the current in rescaled co-ordinates at (a) 
T = 1.23, (b) T = 9.52 and (c) T = 270.9. The last of these is indistinguishable from the steady 
state shape that is found at long times in these rescaled variables. The similarity solution for the 
steady shape in the interior is given by Y = (9X/2) 1 / 3 (dashed line) and is valid away from the 
source and the front regions, which in these rescaled variables requires that T _1 <C X/T <C 1. 




scale, (g/0) 1 / 3 , enters the problem. We thus define a new set of dimensionless variables 
T = t/t* , H = h/h*, X = x/x* and Y = y/y* where t* is an arbitrary timescale and 

x*=y*=( — - — ^ t*. h*=x*ta,n0. (3.5) 
\ q> tan 8 J 

In these non-dimensional variables, the system becomes 

r)fj\ 

V-(HVH)-^), (3.6) 
along with volume conservation in the form 

X„ rY e (X) 

/ H dY dX = f 3 , (3.7) 

where i5 = ^(^tanfl/q) 1 / 3 is essentially the ratio of the two velocity scales in the problem. 
By substituting H = T$(^,?y) with X = T£ and Y = Tij, time can be eliminated from 
this problem completely so that $ is the solution of the two-dimensional problem 

3$ =[«(£ + - 1))] £ + + t])] v , (3.8) 

(with subscripts denoting differentiation) and 

$ d?? d^ = 1. (3.9) 

The system (|3.8ll and l|3.9|l was solved by timestepping the problem in (|3.6|l and l|3.7|l 

using a minor modification of the code described in section 12.31 This was found to be 
a convenient method of solution and also demonstrates that time-dependent solutions 
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converge on the time-independent solution. The results of this calculation are shown in 
figure [3] for a number of different values of 5. 

The importance of the case a = 3 as a transition between qualitatively different flow 
regimes is reminiscent of earlier work on gravity currents. For an axisymmetric gravity 
current, !HupperlHl982l> found that viscous forces dominate inertia at long ti mes for a < 3 
(being insignificant at short times) with the situation reversed for a > 3. lActon et all 
l)200l|) found that a viscous gravity current propagating over a permeable medium spreads 
only a finite distance if a < 3 but spreads indefinitely for a > 3. Despite these similarities, 
the reappearance of a transition at a = 3 here is purely coincidental. 



3.4. a > 3 

In section we observed that for a > 3 a scaling analysis suggests that we should 
observe axisymmetric spreading for T> 1. For such values of a, therefore, we expect to 
recover the axisymmetric solutions given bv lLvle et al\ il2005|) in our numerical simula- 
tions. In particular, for a = 4 we would expect to find that 

X n ,Y max ~ 0.8855T 5 /\ 

where the prefactor here has been determined by repeating the analysis of iLvle et all 
()2005|) . As shown in figure^ this result is indeed obtained from our numerical results. 



4. Experimental results 

We conducted a series of experiments in which a saline solution (dyed red) was injected 
at constant flux (a — 1) into the base of a porous medi um saturated wit h fresh water. 
The details of the experimental setup are as described bv lLvle et all ll2005|) . In summary, 
the experiments were performed in a square-based Perspex tank of internal side length 
61 cm and height 41 cm. The porous medium consisted of a self-supported matrix of 
Glass ballotini (diam eter 3 mm) , which filled the tank to a height of 25 cm. In contrast 
to the experiments of lLvle et all l)2005|) . the Perspex tank was tilted (so that the gravity 
current was propagating on a slope) and the saline solution was injected at the edge of 
the tank, away from the corner because the inherent symmetry is different here to that of 
the axisymmetric case. Video footage of the motion was captured using a CCD camera 
and measurements of the front distance down slope x n as well as the maximum lateral 
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Figure 4. Numerical results for the positions of the current edge X n and Y max as a function 
of time T for a = 4 (solid lines). For T S> 1 these obey the axisymmetric spreading rela- 
ti onship, X„, Ymax ~ 0.8855T 5 / 4 (dashed line), that we expect from the axisymmetric analysis 
of lLvle et cdlfiMh . 



Expt. 


Symbol 


9 (cm s ) 


q (cm 3 s _1 ) 


e n 


r (s) 


x* (m) 


h* (m 


1 


A 


91 


2.14 


9.5 


40.5 


0.112 


0.019 


2 


a 


99 


1.31 


10 


25.2 


0.080 


0.014 


3 


O 


99 


3.04 


18 


11.9 


0.067 


0.022 


4 


• 


99 


4 


18 


13.7 


0.077 


0.025 


5 


■ 


99 


5.78 


18 


16.5 


0.093 


0.030 


6 


★ 


91 


3.86 


5 


196.2 


0.286 


0.025 



Table 2. Parameter values investigated in the six experiments presented here as well as the 
symbol used to represent their results in figure [K] 



extent of the current y m ax were made using the image analysis software ImageJf . The 
details of the six different values of g' , q and 9 investigated are given in table [21 along 
with the relevant values of the typical scales <*, x* and h*. The la tter estimates are based 
on the measurements o f ' 4> = 0-37 and k = 6.8 x 10 -9 m 2 given bv lLvle et all l|2005l) . The 
experimental results of lLvle et al\ l|2005j) are in very good agreement with theory once 
the additional factor of <j> in (|2.5|l is included. We therefore believe these values of <fi and 
k to be correct. 

The experimental results plotted in figure [S] shows that the experimental results are in 
good agreement with the theoretical results produced by solving l|2.8|l . The comparison 
between experimentally observed current profiles and those predicted from theoretical 
solutions of (|2.8J) shown in figure is also favourable — particularly away from the 
source region. Two possible mechanisms may account for the slight discrepancy between 
experiments and theory observed: the drag exerted by the solid substrate on the current 
and the fact that the pore Reynolds number in our experiments is typically 0(5). Such 
a value of the pore Reynolds number suggests that we may be ap proaching th e regime 
where Darcy's law begins to break down, which is around Re = 10 llBearlll988l) . 

| ImageJ is distributed by the National Institutes of Health and may be downloaded from: 
http : //r sb . inf o . nih . gov/ i j / 
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Figure 5. Numerical (solid line) and experimental (points) results for the position of the nose 
of the current, X n , and the maximum horizontal extent of the current, Y max , as functions of 
time. The symbols used to represent each experimental run are given in tabled 
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Figure 6. Comparison between the current profiles (dark) observed in Experiment 3 and 
those predicted by numerical computation (super-imposed white lines). 



5. Geological relevance 

Our experimental and numerical analyses have shown that shortly after the initia- 
tion of a constant flux gr avity current (a = 1 ) it begins to spread axisymmctrically in 
the manner described bv iLvle et al\ l|2005h . However, at times much longer than the 
characteristic time t* given in l|2.6|> . the current loses its axisymmetry and propagates 
predominantly downslope. Since it propagates at constant velocity in this regime, the 
current propagates much further and faster in this case than would be the case if it 
remained axisymmetric. This is potentially problematic in a range of practical applica- 
tions, such as the sequestration of carbon dioxide in which super-critical carbon dioxide 
is pumped into a quifers. Since the dens ity of the liquid carbon dioxide lies in the range 
500 ±150 kgm" 3 llchadwick et oi.ll2005|) . it remains buoyant with respect to the ambient 
water and so will rise up any inclined boundaries. 

The time-scale, t*, over which asymmetric spreading develops is of interest to those 
wishing to predict the course of the released current. While it is difficult to evaluate t* 
in a precise manner because of the uncertainties in the properties of the surrounding 
rock, we can perform some estimates on the b asis of the available data from the Slcipncr 
field l|Bickle et q?Jl2005l Iciiadwick et aJ2005|) . In this Norwegian field, around 10 9 kg of 
liquid CO2 is currently pumped into the local sandstone each year. Presumably due to 
geological complications, this single input flux is observed later to separate into around 
ten independent currents propagating within different horizons of the permeable layer, 
each of which has a volume flux lying in the region 0.002 < q < 0.03 m 3 s _1 . Combined 
with typical measured values for the porosity and permeability of 0.7 < k < 5 x 10~ 12 m 2 
and = 0.31 ± 0.04 as well as the C0 2 viscosity, n — 3.5 ± 0.5 x 10~ 5 Pas ijBickle et all 
2005) we can estimate upper and lower bounds on the value of t*. When 9 = 1°, we find 
that 0.03 <t*< 14.2 years. This suggests that the effects of non-axisymmetric spreading 
may indeed be important in the field. Because of the variety of values of the slope that 
we might expect to encounter in any geological setting, we note also that for 9 ^ 1, 
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t* ~ #- 4 /( 3 -«). For constant pumping rate (a = 1), this gives t* ~ 6~ 2 : i.e. the precise 
value of the timescale over which the current becomes asymmetric depends sensitively 
on 9. This suggests that the different spreading regimes discussed here may be observed 
in the field and may also have practical implications. 

Since injection occurs into confined layers of sediment, estimates for the vertical scale 
of the current, h* , are also important. Interestingly, h* is independent of 8 for 9 « 1 
(measured in radians) and a = 1 so that, with the parameter values given above, we 
find 1.2 < h* < 25 m. This suggests that, near the source, the depth of the sediment 
layer may be similar to that of the current (and so exchange, confined flows may become 
significant). However, we expect that the scaling H <~ T~L 3 valid away from the source 
ensures that the present study will remain valid downstream. 

We are grateful to John Lister for access to his code for a viscous current on a slope 
and to Robert Whittaker for discussions. Mike Bickle, Andy Chadwick, Paul Linden and 
John Lister also provided valuable feedback on an earlier draft of this paper. 
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